#!/usr/bin/perl
# 
# Exercise 12.5
# 
# Perform two BLAST searches with related sequences. Parse the BLAST output of 
# the searches and extract the top 10 hits in the header annotation of each search. 
# Write a program that reports on the differences and similarities between the two 
# searches.
#
# Answer to Exercise 12.5
#
# As usual, there are many ways one could answer this problem.  We'll implement the set
# operations "intersection" and "difference", which will show the sets of hits that
# are the same in the two sets, and that are only contained in one or the other of
# the sets but not both, respectively.

use strict;
use warnings;
use BeginPerlBioinfo;     # see Chapter 6 about this module

# declare and initialize variables
my $beginning_annotation = '';
my $ending_annotation = '';
my %alignments = (  );

# These files were saved from blasts performed at the NCBI web site (www.ncbi.nlm.nih.gov)
# and saved in text format

my $file1 = 'XM_040728blast';
my $file2 = 'XM_040737blast';

#
# Get first array of ids of alignments
#
parse_blast(\$beginning_annotation, \$ending_annotation, \%alignments, $file1);

my @hits1 = get_hit_summary($beginning_annotation);


#
# Get second array of ids of alignments
#
parse_blast(\$beginning_annotation, \$ending_annotation, \%alignments, $file2);

my @hits2 = get_hit_summary($beginning_annotation);

#Just save the top ten hits

@hits1 = @hits1[0 .. 9];
@hits2 = @hits2[0 .. 9];

print "1st array:\n";
print join("\n", @hits1), "\n";
print "\n";
print "2nd array:\n";
print join("\n", @hits2), "\n";
print "\n";

#
# Calculate the set intersection of the two arrays
#

print "Intersection:\n";
my @intersection = intersection(\@hits1, \@hits2);
print join("\n", @intersection), "\n";
print "\n";

#
# Calculate the set difference of the two arrays
#

print "Difference:\n";
my @difference = difference(\@hits1, \@hits2);
print join("\n", @difference), "\n";

exit;

# get_hit_summary -
#  extract the summary list of significant alignments from the beginning annotation
# just report the ID for the alignment
sub get_hit_summary {
	my($beginning_annotation) = @_;

	my(@data) = split(/\n/, $beginning_annotation);
	my(@results) = ();
	my($flag) = 0;

	# gather the alignment summary lines
	foreach (@data) {
		# this follows the alignment summary lines
		if(/^ALIGNMENTS/) {
			last;
		# gather the lines
		}elsif($flag) {
			# ignore the empty lines
			next if /^$/;
			# extract the id (second field, delimited by | character)
			my $id = (split /\|/)[1];
			# save the ID from the line
			push(@results, $id);
		# this precedes the alignment summary lines
		}elsif (/Sequences producing significant alignments/) {
			$flag = 1;
		}
	}

	return sort {$a <=> $b} @results;
}

# set intersection of two arrays
sub intersection {
	my($a1, $a2) = @_;
	# $a1, $a2 are pointers to arrays

	my %seen = ();
	my @results = ();

	foreach (@$a1) { $seen{$_} = 1; }

	foreach (@$a2) {
		push(@results, $_) if defined $seen{$_};
	}

	return @results;
}

# set difference of two arrays
sub difference {
	my($a1, $a2) = @_;

	my %intersection = ();
	my @results = ();


	foreach ( intersection($a1, $a2) ) { $intersection{$_} = 1; }

	foreach ( @$a1, @$a2 ) {
		push(@results, $_) unless defined $intersection{$_};
	}

	return @results;
}
